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Unimolecular evaporation in rotating atomic clusters is investigated using phase space theory 
(PST) and molecular dynamics simulations. The rotational densities of states are calculated in the 
sphere+atom approximation, and analytical expressions are given for a radial interaction potential 
with the form —C/r p . The vibrational densities of states are calculated using Monte Carlo simula- 
tions, and the average radial potential at finite temperature is obtained using a recent extension of 
the multiple range random- walk algorithm. These ideas are tested on simple argon clusters modelled 
with the Lennard- Jones interaction potential, at several total energies and angular momenta of the 
parent cluster. Our results show that PST successfully reproduces the simulation data, not only 
the average KER but its probability distribution, for dissociations from LJ14, for which the product 
cluster can effectively be considered as spherical. Even for dissociations from the nonspherical LJs, 
simulation results remain very close to the predictions of the statistical theory. 



I. INTRODUCTION 

Fragmentation in finite systems offers a convenient way 
to investigate their physical and chemical properties. In 
this respect atomic and molecular clusters have received 
a great deal of attention, and experimental measurements 
of structural 1-7 or electronic 8,9 data have been reported 
using unimolecular dissociation analyses. In particular, 
the relative stability of a cluster is commonly charac- 
terized by its dissociation energies, giving rise to the 
well known "magic numbers" in the mass spectra. 10 Re- 
cently, fragmentation has been employed to probe ther- 
modynamical properties on a more global scale, with 
a focus at phase transitions. Habcrland and cowork- 
ers have used photoabsorption induced fragmentation 
to extract caloric curves in charged sodium clusters ac- 
cross the solid-liquid phase change. 11 They have also 
extended their measurements to probe the liquid-vapor 
phase change. 12 At the same time, Gobet and cowork- 
ers used event-by-event data analyses of multifragmenta- 
tion in Rg (H2) m clusters induced by collisions with a he- 
lium target, 13 showing a small backbending in the caloric 
curve. Brechignac et al. also found some evidences of the 
liquid-gas transition in small strontium clusters from the 
shape of the kinetic energy release distribution subse- 
quent to photoexcitation. 14 

The possible correlations between statistical frag- 
mentation and phase transitions have found theoreti- 
cal supports in cluster physics, 15-18 but also in nuclear 
physics. 19 Decaying nuclei resulting from collisions 20 typ- 
ically show features due to a very large energy deposit, 
where multiple fast fragments are emitted on a short time 
scale. In this case, the main concern is to characterize the 
distribution of the fragments, and the size of the remain- 
ing droplet. In suitable situations, Fisher's formula 21 
gives a correct account of the mass distribution measured 
in experiments. 



Atomic clusters are usually treated much more gently, 
by adding a small amount of excitation energy. Only a 
very few atoms undergo dissociation, and the evapora- 
tive process can take place over long time scales. For ex- 
ample, large weakly bound rare-gas clusters can exhibit 
extremely small rate constants if their excitation energy 
lies not far above the dissociation threshold, because the 
time required for the excitation energy to be located on 
the few dissociative modes rises sharply as the cluster 
size increases. In these systems, one is more interested 
by a complete characterization of the evaporation event 
itself, with a single ejected atom involved. Two observ- 
ables carry most of the useful information, namely the 
dissociation rate and the kinetic energy release (KER) 
distribution. Weerasinghe and Amar (WA) 15 theoreti- 
cally investigated in great details the evaporation pro- 
cess in small argon clusters. Their results show that 
the evaporation rate, and even more the average KER, 
can be used as a probe of the solidlike-liquidlike phase 
change in the parent cluster. 15 To achieve this result, 
they compared various statistical theories of unimolecu- 
lar dissociation to the outcome of molecular dynamics 
(MD) simulations, below the energy range where MD 
becomes prohibitive. One of their conclusions is that 
phase space theory (PST), in the sense of Chesnavitch 
and Bowers, 22 is able to describe accurately the full evap- 
oration statistics in Ar„ clusters, while simpler theories 
such as the Rice-Ramsperger-Kassel (RRK) model 23 or 
the Weisskopf-Engelking formula 24,25 only produce cor- 
rect orders of magnitude. 15 Some important features, in- 
cluding the nonlinear variation of average KER with in- 
creasing excitation energy, are completely absent from 
the predictions of these approximate models. Compara- 
ble methods have been applied by Peslhcrbc and Hase to 
the dissociation in small aluminium clusters, 26 leading to 
similar conclusions. 

Brechignac and coworkers recently reported time-of- 
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flight mass spectrometry measurements of evaporating 
Na+ clusters. 27 A careful interpretation of these results 
necessitated the partitioning of the translational and ro- 
tational kinetic energy released, because only the for- 
mer is actually measured. The possible angular momen- 
tum of the parent cluster is a problem, since rotation 
can strongly alter the evaporation dynamics, hence the 
statistical observables. Up to now, only few theoreti- 
cal works have been devoted to the dynamics of rotating 
clusters. Structural properties and angular momentum 
driven isomerizations were first investigated by Jcllinck 
and Li. 28-30 Using simple statistical theories, Miller and 
Wales further investigated static and evaporation proper- 
ties on the effective rovibrational potential energy surface 
(PES). 31 The influence of angular momentum on clus- 
ter thermodynamics 32 and chaotic dynamics 33 also re- 
ceived some attention. Evaporation in rotating clusters 
had been previously investigated by Stace using simple 
models in the framework of phase space theory. 34 The 
calculations made by this author showed that angular 
momentum tends to increase in small clusters after evap- 
oration (rotational heating), while it tends to decrease 
in large clusters (rotational cooling). These effects have 
been partly observed in the MD simulations performed 
by Weerasinghe and Amar. 15 

As seen from the success of phase space theory to de- 
scribe evaporation in nonrotating clusters, 15 it is highly 
desirable to extend this work to the case of finite an- 
gular momenta. This is the goal of the present paper. 
In the next section, we give the basic PST formalism 
needed to calculate the rotational density of states, in 
the sphere+atom approximation. Exact results are ob- 
tained for a radial interaction potential having the form 
—C/r p , and we provide further details about the numer- 
ical implementation of the method in the more general 
case of a —C/(r — r ) p interaction. The other compu- 
tational ingredients include the estimation of the vibra- 
tional density of states as well as the radial interaction 
potential. We also carry out some MD simulations to be 
used as a benchmark for testing the predictions of PST. 
Application is made in Sec. IV to the evaporation in Ari4 
and Ars, modeled using the common Lennard- Jones po- 
tential. We finally summarize and conclude in Sec. V. 



II. PHASE SPACE THEORY 

In this Section, we work out the main expressions for 
the distribution of kinetic energy released during evapo- 
ration of rotating polyatomic molecules. Conservation of 
angular momentum J is rigorously included in the phase 
space theory 22 ' 35,36 This is particularly important when 
treating rotating systems with prescribed values of J. 
Additionally, PST is built upon the hypothesis of a loose 
transition state, i.e. the products are the transition state. 
In rotating clusters, the centrifugal barrier previously ex- 
plicitely considered by Miller and Wales, 31 is naturally 
accounted for in PST. 



Here we consider a parent cluster characterized by a ro- 
tational angular momentum J and a total rovibrational 
energy E. We denote by J r the rotational angular mo- 
mentum of the subcluster (product) after dissociation. 
Following WA and Jarrold, 37 the probability of finding 
a dissociation event with e tr kinetic energy released is 
given within de tr by: 



i-E—E}. 



P(e tI ,E,J) = R(e tI ,E,J) 
with the differential rate R(e tr , E, J) 



R(£tr, E, J)d,Str , 
(1) 



R(sti, E, J) = Rq 



n ( n J) (E-E ( J) 



e tI )T(e tI ,J) 



(2) 



In the latter equation, Rq is a constant factor that ac- 
counts for channel and rotational degeneracies. 15 E r is 
the rotational energy of the parent cluster. fi^+j and 

Sln^ are the vibrational densities of states (VDOS) at 
angular momentum J of the parent and product clus- 
ters, respectively. T is the rotational density of states 
(RDOS) of the fragments. In these notations, we have 
implicitely assumed that both densities of states of the 
product cluster depend on J, but depend only weakly 
on J r . In the same line of ideas, we consider that the 
clusters are large enough so that the energy difference 
Eq'^ between the potential energy minima of the parent 
and product clusters can be taken at the same value of J 
for both clusters. The knowledge of the differential rate 
function R readily leads to the average kinetic energy 
released: 



<£t r > = 



j-E—E^ 

J min 



e t iP(e t i,E, J)de tT . 



(3) 



The calculation of P(s tT ,E, J) and (etr) requires one to 
compute both il and T for the product cluster, but nei- 
ther the constant Rq nor the VDOS for the parent clus- 
ter. The absolute rate constant obtained from integrating 
R(str, E, J) with respect to e tr thus requires a more sub- 
stantial effort than any data related to the KER. In the 
remainder of this section, we focus on the rotational den- 
sity of states. The vibrational quantities are relatively 
easy to compute, and they will be dealt with later. 

The calculation of the T(e tr ,J) function is directly 
linked to the energetics and angular momentum con- 
straints during dissociation. Our main assumption will 
be to treat the evaporative system LJ n+1 -^LJ„+LJ as 
well represented by a sphere+atom model. Then the 
product cluster LJ„ has a unique rotational constant B. 
In this case, Chesnavitch and Bowers have shown that 
the rotational density of states can be calculated as 22 



r(etr, J) 



r(e*, J r )dJ r dL, 



(4) 



where e* stands for the upper limit of the rotational 
energy, L for the orbital angular momentum of the 
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products. The integration is carried out in the (J r ,L) 
plane with the boundaries S discussed below. In the 
sphere+atom case, T(e*,J r ) simply equals 2J r , 37 hence 
the problem is reduced to finding the expressions for the 
boundaries S. We first consider the case of a radial dis- 
sociation potential V(r) given by V(r) = —C/r p , with p 
greater than 2. 

The first boundary on S is given by the constraint on 
the kinetic energy of the dissociating atom. The centrifu- 
gal barrier that must be overcome is located at i 



such that Vl(t) 
yields 



V(r) + L 2 /2[ir 2 is maximum. This 



E t =i Vb-2)/ Ap) 



with the notation 



A,, 



_ C 2/(p-2) 



P 



h 2 



p/(p-2) 



(5) 



(6) 



For the atom to actually dissociate, its kinetic energy 
must be positive at the barrier, which is expressed as 



SJ r 2 + L 2 ^- 2 )/A p < £tr . 



(7) 



The second boundary on S comes from the conservation 
of angular momentum, J = J r + L, or 



\J r — L\ < J < J r + L. 



(8) 



The conditions (7) and (8) define lower and upper bounds 
for J r at each value of L, denoted as J™ m (L) and 
J™(L), respectively. Integration over the contour S 
can be formally carried out: 



T(etr, J) 



-L [<J " 



where we have introduced the lower and upper bounds for 
the integration on L, namely L m i n and £ max - Let now 
C be the set of (J r ,L) points that fulfill the equation 
£ tr = L 2 P/(P- 2 )/A p + BJ 2 . Let also J r * and L* be the 
intersection points of C with the abscissa and ordinates 
axes, respectively. We find that 



J* = (str/B) 1 / 2 and L* = (A p e tr ) 



(p-2)/2p 



(10) 



Four different cases must be treated separately, depend- 
ing on whether the values of J* and L* are smaller or 
larger than the initial angular momentum J. These four 
cases are depicted in Fig. 1. They correspond to the fol- 
lowing conditions: 

(a) J < J*; J < L* , 

(b) j;<j<l\ 

(c) L*<J<J*, 

(d) J r * <J;L*< J. 

In cases (a) and (c), integration starts at £ m j n = 0. In 
cases (b) and (d), L m i n is determined by the intersection 
of C with L = J — J r . We denote this point by L\\ 



L J J T 1 

L 2p/(p-2)/ Ap . 



BJ 2 



(11) 



In (a) and (b), L max = 
intersection of C with L 



Li is obtained at the unique 

-- J + J r : 



L = J+ J r , 

L 2 P /( P -2) /Ap 



BJ 2 



£tr 



(12) 



Finally, the upper bound L max = L 2 is given in cases 
(c) and (d) by the intersection of C with L = J — J r , 
equation (11) above. Thus in (d) the two extremal values 
are solutions of the same equation. 

At any value of L in the range L m [ n < L < L max , 
the lower and upper values of J r (L) also depend on the 
conditions (a-d). For example, in case (a) one must dis- 
tinguish between three subcases, namely < L < L\, 
Li < L < J, and J < J < L 2 , where L\ denotes 
the intersection of C with L = J r — J. In the range 
< L < J, J™ in is equal to J - L. In the range 
J < L < L 2 , Jf n {L) = L - J. The upper bound 

jmax ig giyen by jmax = J + ^ for < i < Li , and 

by J r max (L) = [e tr - L 2 p/(p- 2 ) /A p ]/B for L x < L < L 2 . 

After some algebra, integration over the boundary S 
leads to the total rotational density of states: 



r(e tr ,J) = (L 2 — Li)(e tI /B — J 2 



A P B 3p-2 



— L, 



+J(I 2 + Ll)-{Ll-L\)/Z. (13) 

We will not discuss the three other cases (b-d) in details, 
and we only provide below the final results. In cases (a) 
and (c), the RDOS is given by Eq. (13) above. In cases 
(b) and (d), it is expressed by 

r( £tr ,J) = (L 2 — L 1 )(e tI /B — J 2 ) 



1 „ _ O / 3p-2 3p-2 
1 I' ~ I J T T ~T 



A P B 3p-2 



+J(L 2 -L\)-(Ll 



L?)/3, (14) 



which only differs from Eq. (13) by the quantity 2JL 2 . 

We have not yet discussed the lower bound of integra- 
tion on e t r in Eqn. (1) and (3), denoted as £™ m . This limit 
occurs when the curve C is tangent to the line L = J — J r . 
This condition can be cast into an equation in J r only: 



p-2 



BApJ r (J — J r )'- 



1. 



(15) 



which can be easily shown to have a unique solution in 
the range < J r < J. The value of e" r lln follows from 
substitution in (C) and L = J — J r . In the case of a 
van der Waals dispersion interaction, V(r) = —C/r 6 , the 
appropriate values for the rotational density of states are 
given exactly for p — 6 by 



r(e tr ,J) = (L 2 — Li)(e tr /B - J 2 ) 

+J{Ll±L\)-{Ll-L\)/Z 



^2 _ ^1 

4ArB 



(16) 
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FIG. 1: Schematic representation of the (L, J r ) integration plane, denoted by S in the text, (a) J < L* and J < J*; (b) 
J* < J < L*; (c) L* < J < J*; and (d) J > L* and J > J*. The outer boundary is defined by the curve C (see text), the 
inner boundary is the circle with radius J centered at (0,0). 
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where the plus sign stands in cases (a) and (c), and the 
minus sign stands in (b) and (d). In the case p — 6, an 
anaytical expression for £™ ln can also be found 



BJ 2 + -5 2 A 6 J+ ^5 3 A 2 



2 

27" 



2 

27" 



^B 3 A 2 + lB 2 k*J 



6J 
BA* 



1/2 



which is given in the low J regime by £™ m = J 3 /A 6 , up 
to the fourth order in J. The numerical implementation 
of the above formulas is straightforward. At fixed total 
energy E and angular momentum J of the parent cluster, 
one must first calculate the lower and upper limits e™ m 

and e™ ax = E — E^ J \ For each value of e tr in this range, 
equations (11) and (12) must be solved numerically to 
give Li and L 2 , and the rotational density of states is 
calculated using Eqn. (13) and (14) above. 

All the previous formalism has been derived by assum- 
ing an interaction potential with the form V(r) = —C/r p . 
It turns out that this expression does not give a very good 
account of the finite extension of the cluster, and that a 
better representation of the atom-cluster interaction is 
provided by V(r) = —C/(r — r ) p , with r > 0. In this 
case, and more generally for an arbitrary form of V(r), 
the computation of the rotational density of states must 
be carried out numerically. Firstly, for a series of L, the 
location r*(L) of the centrifugal barrier is obtained by 
solving 



dV 
dr 



L 2 



*\3 ' 



(18) 



Actually this dependence acts by two ways. 32 Firstly, the 
centrifugal effects perturb the potential energy surface 
into an effective, rovibrational surface. 28 ' 31 Secondlya, 
the conservation of the vector J adds an extra geomet- 
rical weight in the configurational density of states or in 
the pa rtition function. This weight is given explicitely 
by l/Vdetl, where I is the inertia tensor at the current 
configuration. 32,39 The angular momentum is naturally 
conserved in constant energy molecular dynamics simu- 
lations, as well as in constant temperature Nose-Hoover 
schemes at J = 0, but not in conventional Monte Carlo 
simulations. Therefore, some differences between the MD 
and MC procedures may arise when this weight may span 
several orders of magnitude, as in the case of Ar^ , which 
is linear in its ground state geometry. 40 

To calculate fl^(E) for several values of J, we have 
used the Monte Carlo method proposed in Ref. 32, fur- 
ther improved with the parallel tempering accelerating 
scheme. 41 The multiple histogram method 42 was then 
used to estimate the configurational densities of states. 
In turn, the total vibrational densities of states were ob- 
tained from the configurational densities by a simple con- 
volution product. 

The present results were checked by performing ad- 
ditional molecular dynamics simulations. The histogram 
analysis 38 showed a good agreement between the MC and 
the MD calculation. We have also checked the physical 
relevance of the calculation by computing other related 
thermodynamical observables, such as the canonical heat 
capacity. At low values of J, the melting temperature 
was seen to roughly decrease with J as J 2 , in agreement 
with previous works. 32 



One then deduces the height (L) of the barrier: 

L 2 



e\L) 



V(r*)- 



2u{r*) 2 

At a given e tr , the integration boundaries S become 



(19) 



e\L) + BJ 2 <e tr , 
\J r - L\ < J < J r + L 



(20) 



and the limits L\ and Li are still given by equations 
(2) and (3) after replacing L 2 p/(p- 2 )/A by e*(L). While 
(J™ 11 ) 2 is still equal to (J-L) 2 , (J™ ax ) 2 is now obtained 
from [e tr — e^]/B in the range L\ < L < L 2 - The integra- 
tion of r(e tr , J) must be done numerically, after estimat- 



ing e™ m from the tangency condition of £tr 
with L = J r — J. 



£ T + BJ l r 



III. COMPUTATIONAL PROCEDURE 

A. Vibrational density of states 

The vibrational densities of states Q( j \E) depend im- 
plicitely on the total angular momentum J of the cluster. 



B. Radial potential 

In a first approach, the dissociation potential V(r) felt 
by an atom leaving the n-atom Lennard-Jones cluster 
can be approximated by its asymptotic form —C^/r 6 . 
At very large distances r, the parameter is given 

by Ansa 6 LJ units. However, at intermediate distances, 
where the centrifugal barrier is likely to be located, the fi- 
nite extent of the cluster induces significant deviations of 
the average potential. A simple approach to this problem 
is to consider a continuous homogeneous distribution of 
Lennard-Jones centers inside a sphere of radius R oc n 1 / 3 . 
For large sizes, this leads to the Gspann-Vollmar poten- 
tial V n A3 



V n (r) = Cia 



r 6 + 21r 4 r 2 /5 + 3r 2 r 4 + rg/3 



(r 2 - r 2 f 



C 6 



(r 2 — r 2 ) 3 ' 

where C12, C§ and tq are size-dependent and given by 



(21) 



C12 = Ansa ; 
C 6 = Ansa 6 ; 



(22) 



r = {3/pnp) 1 ^ 3 - 1]. 



(23) 



0.0 



In the above equation, p is the atomic density in the 
solid state. The Gspann-Vollmar potential was built in a 
similar way as the Girifalco potential describing the in- 
teraction between Ceo molecules. 44 It is not appropriate 
for medium-size, nonspherical clusters, or for interme- 
diate distances r, where the continuous approximation 
would break down. In addition, because the cluster is 
thermalized at a finite temperature, the atomic fluctua- 
tions may induce some changes in the average potential 
felt by the tagged distant atom. Using constraint dy- 
namics, Weerasinghe and Amar 15 showed that the sim- 
ple —Cs/r e form was not fully appropriate to describe 
the atom-cluster interaction, and that a much better fit 
was obtained using the — C§j{r — ro) 6 form. The Gspann 
and Vollmar results also suggest that a —Cq/(t 2 — Tq) 3 
form could also be used. We have carried out some con- 
strained Monte Carlo simulations at finite temperature, 
by keeping the external atom at a fixed distance r from 
the cluster center of mass. The temperature effects were 
not investigated in the paper by Weerasinghe and Amar, 
and we have chosen to perform the calculation at low 
(0.01 LJ units) and high T, namely T = 0.2 for LJs and 
T = 0.3 for LJ14. These values are close to the melting 
points of the two clusters, above which evaporation takes 
place in a sub-nanosecond time scale. 

As an alternative to the MC simulations with con- 
straints, we have calculated the finite temperature dis- 
sociation potential using the recently proposed multiple 
range random walk algorithm by Wang and Landau. 45 
This method has been straightforwardly extended to 
the computation of effective potentials and potentials of 
mean forces, 46 and consists of performing a Metropolis 
Monte Carlo simulation using the following acceptance 
rule: 
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acc(R id — » Rnc 
|^exp{-/?[L(R n 



= 

)-L(R old )]}| , (24) 



where (3 = 1/ksT, R id and R n0 w are two successive 
points in the configuration space, r id and r new the cor- 
responding atom-cluster distances, respectively. g(r) is a 
weight function, initially set to 1 in the entire range of 
r, which evolves dynamically along the MC simulation 
by the operation g(r) — > / x g(r) after the distance r 
has been visited. The constant factor / is initially set to 
2, and gradually decreases to 1 after a given number of 
Monte Carlo steps. After several iterations of this pro- 
cess, the function L(r) = —(3~ 1 lng(r) converges to the 
potential of mean force (PMF) W(r): 4e 



T(r) -» W(r) = -fi-Hnpir), 



(25) 



where p(r) is the probability distribution of finding the 
atom at the distance r from the cluster, given by the 
canonical average p(r a ) = (S[r — r(R)]). 

Once the PMF is known, it can be subsequently used in 
a biased multicanonical simulation to sample the entire 



FIG. 2: Atom-cluster radial potential in the reaction 
LJn+i ^LJ„+LJ. The symbols are from constrained MC sim- 
ulations, the solid and long dashed lines are the result of the 
Wang-Landau (WL) multicanonical reweighting scheme. The 
asymptotic law —Cq/t , with C% = An LJ units, is also drawn 
as a dashed line. The data are plotted for two temperatures 
in each case, (a) n = 7; (b) n = 13. 



range of distances in a uniform way. For this we replace 
the potential V by V + W. The average potential V(r) 
felt by the atom at distance r from the cluster center of 
mass is given by the usual reweighting formulas 47 The 
Wang-Landau scheme allows one to compute the poten- 
tial V over a continuous range of r, instead of only a 
small set when using constraint dynamics. 

We have represented in Fig. 2 the effective potential 
calculated from constrained MC simulations and from 
multicanonical simulations for the two sizes n+1 = 8 and 
n + 1 = 14, at low and high temperatures. In each case, 
we observe a very good agreement between the two meth- 
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ods, suggesting that the Wang-Landau/multicanonical 
scheme can yield accurate average potentials over con- 
tinuous range of distances. For the two clusters, the in- 
teraction between the external atom and the cluster is 
stronger at higher temperature. This can be explained 
by the larger "volume" of the cluster in its liquidlikc 
state with respect to its solidlike low temperature value. 
Hence the apparent extent of the cluster is larger, and 
the potential is larger in modulus. As can be seen in 
the two pannels of Fig. 2, the asymptotic — C 6 /r 6 form 
is not appropriate at intermediate distances when tak- 
ing Cq = Ansa 6 . As noted by Weerasinghe and Amar, 15 
setting this constraint on C§ free does not improve the 
behavior of V(r) much. A better fit is obtained with the 
expression V(r) = —Ce/(r — ro) 6 . The values of Cq and 
r as a function of cluster size and temperature are given 
in Table I. 

TABLE I: Fitting parameters Cq and ro of the average atom- 
cluster LJ potential V(r) — ~Ce/(r — ro) 6 , and average rota- 
tional constant B, for LJ„ clusters at different temperatures. 



Cluster 


Temperature 


c 6 


ro 


B 


size n 


(e/fcs) 


(sa 6 ) 


(*) 


(ma 2 ) 


7 


0.01 


8.070 


0.809 


0.150 


7 


0.20 


8.728 


0.802 


0.135 


13 


0.01 


19.149 


0.733 


0.0534 


13 


0.30 


62.353 


0.499 


0.0444 



Temperature effects are weak on the smaller cluster, as 
both Cq and r remains nearly constant. This is probably 
due to the fact that, upon melting, the non spherical LJ7 
cluster does not really enlarge, but instead visits other 
(still non spherical) isomers. On the other hand, a sig- 
nificant change is seen on the parameters for the much 
more spherical n = 13. At large temperature, the val- 
ues we get are found to yield a very similar average po- 
tential than the one found by Weerasinghe and Amar. 15 
However, because these authors employed constant en- 
ergy MD simulations and because they did not provide 
the total energy used for this cluster, we cannot reliably 
compare our results with theirs. 

We now have all the ingredients required for the PST 
calculations. In order to assess or question the quality 
of the statistical theory, we need to carry out simula- 
tions of the actual evaporation process at finite angular 
momentum. 

C. Molecular dynamics simulation of the 
evaporation dynamics 

For each size and for each value of angular momentum, 
we have considered a set of 5000 molecular dynamics tra- 
jectories. The average kinetic energy release has been 
analysed after each trajectory ending into an evaporation 
event. The two contributions of the KER were evaluated, 
namely the rotational part of the product LJ„ subcluster, 



and the translational contribution of the external atom 
undergoing evaporation. The instant of evaporation was 
considered as the last time at which the radial velocity 
of the atom was negative. 15 By varying the total angular 
momentum between J = and J w 5 LJ units, 48 we have 
performed a systematic study of the effects of rotation on 
the evaporation process. 



IV. RESULTS AND DISCUSSION 

In this paper, we will focus on the energetics of the 
unimolecular process, rather than on the absolute evap- 
oration constant. This choice is mainly guided from 
previous studies where the relationship between phase 
transition in the product clusters and the evaporation 
statistics was most clearly evidenced on the kinetic en- 
ergy released. 15-17 Two different unimolecular reactions 
have been considered, involving the nearly spherical LJ13 
product and the nonsphcrical LJ7 system. This lat- 
ter cluster has an ellipsoidal symmetry with inertia mo- 
menta in the ratio (0.64,0.64,1) at T = 0. The two 
reactions studied here are thus LJ14 — >LJi3+LJ and 
LJ 8 -^LJ 7 +LJ. 



A. LJ14 — >LJi3+LJ 

LJ13 is a magic cluster and shows enhanced melting 
point and latent heat of melting with respect to its im- 
mediate neighbors LJi 2 and LJ14. In addition to its 
spherical shape, this cluster provides a good candidate 
for investigating the melting transition as seen from its 
evaporation observables. 

To apply the PST formalism, we first have to calcu- 
late the rotational density of states T(s tVl J) for different 
values of the angular momentum J. As explained in Sec- 
tion II, r(e t r, J) is sensitive to the radial potential V(r) 
felt by the dissociating atom. It is also sensitive to the 
rotational constant B of the product LJ13 cluster. In 
Fig. 3, r(e t r, J) is plotted for J=0, 2, and 4 using the ra- 
dial potential —C§/(r — r ) 6 and the rotational constant 
calculated at moderate temperature, T = 0.3, close to the 
melting point. The parameters are given in Table I. For 
comparison, we have also plotted the RDOS calculated 
using the simpler radial potential — Cg/r 6 , with = An 
and the rotational constant at T = 0. The latter poten- 
tial does not fully reproduce the simulated potential, but 
the asymptotic limit (large r) is known to be exact. Be- 
cause the discrepancy between the two radial potentials 
is quite large, the significant differences between the pre- 
dictions of PST in the rotational densities show that the 
centrifugal barriers are indeed located within the range 
of distances plotted in Fig. 2. The PST calculation using 
the analytical results at r = notably underestimates 
the RDOS, by about 15% for any J and e tr . The ef- 
fects of a finite angular momentum in the parent cluster 
are strong. The value of etr at which T sharply increases, 
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FIG. 3: Rotational density of states T(etr, J) as a function 
of £ t r for 3 values of J in the unimolecular dissociation of 
LJ14. The curves plotted are the predictions of PST using 
the simulated radial potential at T = 0.3 (thick lines) or the 
simpler —Ce/r 6 potential with Ce = 52 (thin lines). 



FIG. 4: Average translational ((etrans), squares), rotational 
((e r ot), diamonds) and total ((etr), circles) kinetic energies re- 
leased as a function of J in the dissociation of LJ14 at £=-26. 
All results are from MD simulations. Full symbols correspond 
to E = —26; empty symbols are for E = —29. 



previously denoted as £™ ln , clearly changes with J. More- 
over, the slope of the RDOS also increases with J. As a 
consequence, there is an order of magnitude increase at 
£tr = 3 between J = and J = 4. This latter value can 
be considered as still moderate for LJ14, as there is only a 
small variation in the potential energy surface, and in the 
related properties such as the heat capacity. 31,32 In par- 
ticular, the cluster structure is only slightly perturbed, 
and spontaneous isomerization is not expected to take 
place below about 15 LJ units. 29 

In Fig. 4 we have represented the results of the molec- 
ular dynamics simulations to be used as benchmarks for 
the present theoretical analyses. The rotational, transla- 
tional and total kinetic energy released are plotted as a 
function of J for two different total energies of the parent 
cluster. The effect of J on the energetics of the dissocia- 
tion reaction is mainly governed by the evolution of the 
rotational contribution of the KER. Around J = 4, this 
contribution becomes larger than the translational en- 
ergy. This latter contribution appears almost constant 
up to J ~ 4 and slighty increases at higher J. In the 
range of energies considered here, the effect of internal 
energy is weak. We observe a steady increase in the en- 
ergies released during evaporation as both angular mo- 
mentum and internal energy increase, as intuition would 
suggest. 

In the form detailed above, phase space theory only 
gives us access to the total (translational+rotational) ki- 
netic energy released. We have plotted in Fig. 5 the 
variations of the average KER (e tr ) calculated from MD 
simulations as a function of J, and compared them to 
the predictions of PST in the following three approxima- 
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FIG. 5: Average total kinetic energy release in the dissocia- 
tion of LJ14 as a function of angular momentum J at £=-26. 
Comparison between MD simulation and phase space theory 
using different radial potentials. 



tions. The radial potential was either taken with ro = 
or with finite r , and in the latter case, two tempera- 
tures were taken, corresponding to either the rigid case 
(T = 0.01) or to the liquid case (T = 0.30). The values of 
the parameters ro, Cq and B implicitcly depend on these 
approximations. As can be seen from Fig. 5, all three ap- 
proximations perform quite well. Looking more closely, 
we notice that the PST calculation at low temperature 
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overestimates (etr), and that the approximation tq = 
leads to a slightly diverging KER above J ~ 3. Actually, 
only the calculation at finite ro and temperature close 
to the melting point remains quantitatively close to the 
simulation data in the entire range of E. This is not so 
surprising, because the simulation takes place at rather 
large internal energies, where the cluster is in a liquid- 
like state. At low temperature, the average rotational 
constant is larger (see Table I). Therefore the rotational 
contribution to e tr is overestimated, which explains the 
relatively high values of (e tr ) in Fig. 5. Another possible 
cause could be the less attractive interaction potential 
at this low temperature (see the lower panel of Fig. 2), 
resulting in higher centrifugal barriers, hence a further 
shift of £tr to a larger value. 

The variations of the average KER as a function of in- 
ternal energy are displayed in Fig. 6 for the same three 
values of angular momentum, J = 0, J = 2, and J = 4 
LJ units. The results of the PST calculations are repre- 
sented in the energy range where an inflection occurs due 
to the melting phase change in LJ13. 15 The excitation en- 
ergy where this inflection takes place show a weak depen- 
dence over the radial potential used, as well as a nearly 
constant value with increasing J. As angular momentum 
increases, the average KER also increases, following the 
expected scaling law (£ tI )(J > 0) w (str)(J = 0) + aJ 2 . 
Looking now at the differences between the PST calcu- 
lations, we notice that the use of the — C/r 6 radial po- 
tential overestimates the average energy released, more 
and more as J increases. This quantitative difference can 
be explained from the differences in the RDOS, as seen 
in Fig. 3. For excitation energies E/n <~ 1, (e tr ) is less 
than about 1. In this range, the energy shift e™ ln plays 
a crucial role, and its overestimation with the ro = ap- 
proximation is consistent with the larger average kinetic 
energy released in Fig. 6. 

Because Fig. 5 does not confidently discriminates be- 
tween the various approximations used in our application 
of PST, we turn to the probability distribution of e tr at 
given total energy and angular momentum of the par- 
ent cluster, which carries more precise information. In 
Fig. 7 we have represented these distributions obtained 
from MD simulations and from PST, with the two radial 
potentials cither with tq — and Cg = 4n or with the 
potential extracted from simulations close to the melting 
point. The results are given for J = and J = 5 LJ 
units. As can be observed on the upper panel of this 
figure, the agreement between MD and PST is excellent 
at J = 0, and the two PST calculations give very similar 
data. The rotational density of states roughly shows lin- 
ear variations upon increasing £ tr in Fig. 3. We also show 
in Fig. 5 the probability densities computed using the ex- 
plicit linear approximation r(e tr , J) = a(e t r — £ ti m )- Ac- 
tually the slope constant a does not play any role in the 
KER distribution. This linear Weisskopf-like behavior 2 
leads to slight shift of the distribution toward higher en- 
ergies. 

A nonzero angular momentum appears as a more strin- 
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FIG. 6: Average total kinetic energy release in the dissocia- 
tion of LJ14 as a function of E/n for 3 values of the initial 
angular momentum J, from the predictions of PST using the 
simulated radial potential at T = 0.3 (thick lines) or the sim- 
pler — C(i/r e> potential with Ce = 52 (thin lines). 
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FIG. 7: Probability distributions of the total kinetic energy 
release etr at E — —26 for two values of J, in the dissociation 
of LJ14: (a) J=0; (b) J=5. The dashed lines refer to the 
linear approximation r(e tr , J) oc Etr — e™ m (J). 



gent test for the statistical theory, as we see on the lower 
panel of Fig. 7 that the agreement between MD and PST 
is significantly better with the calculation performed us- 
ing the radial potential corresponding to T = 0.3. This is 
consistent with the better agreement previously observed 
in Fig. 4. In particular, the threshold value e™ ln of e tr at 
which the probability suddenly rises (near e tr <~ 0.5 LJ 
unit) is well reproduced by PST at T — 0.3, but is not 
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in the approximation ro = 0. Evaporations from rotat- 
ing clusters are characterized by a nonzero value of this 
threshold, which is due to the extra excitation energy 
needed for the system to overcome the larger centrifugal 
barrier. The influence of the shape of the rotational den- 
sity of states on the probability distribution of e tr is also 
investigated using the linear approximation. The differ- 
ence with previous PST calculations is more significant, 
which indicates the strong influence of the shape of the 
RDOS in the vicinity of e™ in . 
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FIG. 8: Average angular momentum of the product subcluster 
as a function of the initial angular momentum J of the parent 
LJ14, from MD simulations at two total energies. 

Finally we have represented in Fig. 8 the average final 
angular momentum J r of the product cluster as a func- 
tion of the initial J, at two total internal energies, from 
molecular dynamics simulations. On this figure the line 
J r — J = is also drawn. Depending on the initial angu- 
lar momentum, the product cluster can gain or lose some 
of its rotational velocity. These rotational cooling and ro- 
tational heating effects occur for J < Jo = 2.3 LJ units 
and J > J , respectively. The threshold value Jo weakly 
depends on the total internal energy in the rather limited 
range investigated here. 34 This can be understood using 
the following simple arguments. At low J, evaporating 
an atom induces an orbital momentum L, which is nearly 
balanced by the angular momentum J r . Hence angular 
momentum increases for initially small values of J. On 
the other hand, rapidly rotating clusters lose a part of 
their angular velocity by emitting one atom, and J tends 
to decrease upon evaporation. 
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We now turn to the evaporation statistics in the 
smaller LJg cluster. The product cluster LJ7 is nonsphcr- 



FIG. 9: Average total kinetic energy release in the dissocia- 
tion of LJg as a function of angular momentum J at £=-10.21. 
Comparison between MD simulation and phase space theory 
using different radial potentials. 



ical not only in its lowest energy structure, but also in 
any of his three other stable isomers. The average ki- 
netic energy released during dissociation is represented 
in Fig. 9 versus the total angular momentum of the par- 
ent cluster J. The values plotted are the results of MD 
simulations as well as the predictions of PST under var- 
ious approximations concerning the radial potential. As 
in the previous paragraph, we have considered the simple 
—Ce/r e case with Cq = An, and two more realistic poten- 
tials extracted from Monte Carlo simulations at T = 0.01 
and T = 0.2, respectively. The latter value is close to the 
isomerization point in this system. The effective rota- 
tional constant was taken as the average over the different 
instantaneous values at the corresponding temperatures. 
They are also given in Table I. As in Fig. 5, the gen- 
eral agreement between MD and PST is good, but we 
notice that the discrepancy is larger for the PST calcula- 
tion with ro — 0. The effect of temperature on the radial 
potential is weak in this case. This may be partly due 
to the lower melting point of this system, but should be 
correlated with the similar radial potentials felt by the 
leaving atom, as represented in Fig. 2. 

The probability distribution of kinetic energy released 
is reported in Fig. 10 at the total energy E = —10.21 
LJ units, and at zero or nonzero angular momentum of 
the parent cluster. As was the case for the bigger clus- 
ter, PST reproduces very accurately the results of MD 
simulations at J = 0, and the two calculations with the 
different radial potentials yield essentially similar data. 
In contrast, the distributions at J = 3 differs somewhat 
from the simulation results. In particular, the general 
shape predicted by PST is too sharp with respect to 
MD, and the value of e tr where the probability starts 
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FIG. 10: Probability distributions of the total kinetic energy 
release e tr at E = —10.21 for two values of J, in the dissoci- 
ation of LJ 8 : (a) J=0; (b) J=3. 

to increase is too high by about 50%. This error further 
increases when using the alternative radial potential cor- 
responding to r = 0. The discrepancies observed here 
should be mainly due to the erroneous assumption of a 
spherical product. 

Nevertheless, the global behaviors of both the KER 
and its probability distribution remain correctly repro- 
duced by phase space theory, indicating that this sta- 
tistical approach captures all the important physical 
and chemical ingredients of unimolecular dissociation in 
weakly bound systems, especially the conservation of an- 
gular momentum. 

V. CONCLUSION 

Unimolecular decay in large atomic systems is eas- 
ily treated using simple theories such as the RRKM or 
Weisskopf-Engclking statistical approaches. Phase space 
theory not only includes possible anharmonic effects, 
but also the rigorous constraints on angular momentum, 
through general expressions for the differential rates of 
dissociation. Weerasinghe and Amar 15 were the first to 
show clearly that PST is qualitatively and quantitatively 
accurate in predicting rate constants and energetic distri- 
butions in the evaporation of nonrotating atomic clusters. 
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Building upon their seminal paper, we have extended 
their work to the more general case of a finite angular 
momentum in the parent cluster. For this we calculated 
exactly the rotational density of states in the case of an 
interaction between the product cluster and the dissoci- 
ating atom having the form —C/r p . The implementation 
to other forms has also been given. The anharmonic vi- 
brational densities of states were calculated using Monte 
Carlo simulations on the effective rovibrational energy 
surface, 32 and the radial potential was calculated using 
an extension of the recent Wang-Landau algorithm. 45 ' 46 

We have tested the applicability of PST to the case of 
unimolecular evaporation in the LJ14 and LJg clusters. 
These two different cases allow us to question the hypoth- 
esis of a sphere+atom collision underlying the statistical 
formalism. We have shown that PST was quantitatively 
accurate in predicting the distribution and average value 
of the kinetic energy released during dissociation, espe- 
cially after considering the radial potential calculated at 
temperatures close to the melting point, where dissocia- 
tion actually occurs on the time scale of MD. Taking the 
simple form — C/r 6 introduces some extra errors, in par- 
ticular at large energies and angular momenta. Wc have 
also seen that dissociation in LJ 8 was less well described 
using PST in the sphere+atom assumption. Beyond this 
approximation, one could generalize the present formal- 
ism to the case of ellipsoid or even triaxial shapes. Ex- 
tension to molecular systems is also possible, provided 
that the internal degrees of freedom of the dissociating 
molecule are correctly accounted for. 

To bridge the gap between the results obtained in the 
present work and the experimental concerns, one must 
extract the separate translational and rotational contri- 
butions to the total kinetic energy released during evap- 
oration. This separation was achieved previously in the 
case of nonrotating parent clusters. 49 As a first step, it 
would be useful to extend the present effort to the char- 
acterization of the angular momenta distribution after 
dissociation of an initially rotating system. The starting 
distribution of J could there be either a delta function 
(as in the present work) or a thermalized Boltzmann dis- 
tribution P(J) oc J 2 cxp(-BJ 2 /k B T). Work along these 
lines is presently in progress. 
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